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SUMMARY 


This  report  describes  flow  and  temperature  distributions  In  the  Jet 
Fuel  Thermal  Oxidation  Tester  (JFTOT),  a  device  that  has  been  used  In 
thermal  stability  research  of  aircraft  fuels.  A  three-dimensional 
numerical  analysis  Is  presented  In  this  report  using  the  K1VA  code.  This 
work  represents  a  preliminary  attempt  to  "calibrate"  a  Computational  Fluid 
Dynamics  and  Chemistry  (CFDC)  model  using  existing  data  from  a  fiber 
optics  modified  JFTOT. 

Because  of  the  limited  amount  of  data  on  JFTOT  experiments,  the 
deposition  model  used  In  the  study  Is  the  global  type  Arrhenius  equation. 
Upon  the  finding  of  the  parameters  relevant  to  mass  transfer,  heat  trans¬ 
fer,  and  fluid  dynamics  by  future  experiments,  those  relevant  parameters 
can  be  combined  altogether  to  define  thermal  Instability  of  jet  fuel  with 
a  higher  degree  of  accuracy. 
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SECTION  I 


INTRODUCTION 

In  order  to  improve  understanding  of  the  thermal  stability  of  aircraft 
fuels  that  have  been  tested  in  the  Fiber  Optic  Modified  Jet  Fuel  Thermal 
Oxidation  Tester  (FOM-JFTOT) ^ ,  a  three-dimensional  analysis  on  the 
FOM-JFTOT  was  performed.  This  level  of  flow  detail  was  needed  because  of 
the  FOM-JFTOT's  unique  configuration  with  the  exit  tube  rotated  90°  with 
respect  to  the  flow  inlet. 

There  are  perhaps  more  computational  fluid  dynamics  (CFD)  codes^ 
that  have  the  capability  for  three-dimensional  analysis.  All  CFD  codes 
use  one  of  several  methods  to  solve  the  governing  Navier-Stokes  equa¬ 
tions.  Each  of  these  has  its  own  merits  and  shortcomings.  For  this 
effort,  the  three-dimensional,  time-dependent,  finite  difference  KIVA^~® 
code  was  selected  and  modified  to  analyze  temperature  and  flow  distribu¬ 
tions  of  JP-5  fuel  in  the  FOM-JFTOT  owing  to  ease  at  which  this  code  can 
be  used  to  model  different  fluids  flowing  in  arbitrary  geometries. 

Sections  II  and  III  discuss  background  and  general  deposition  models 
shown  in  literature.  Section  II  describes  our  JFTOT  model  and  results 
from  the  study.  Section  V  contains  conclusions  and  recommendations  for 
future  studies. 
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SECTION  II 


BACKGROUND 

Thermal  stability  of  hydrocarbon  fuels^'^  is  expected  to  become 
increasingly  important  in  the  future.  Especially,  coal-derived  jet  fuels 
could  be  crucial  to  our  national  security  if  the  supply  of  imported 
petroleum  is  disrupted.  Speaking  at  the  Spring  National  Meeting,  John  R. 
Rindt  of  the  University  of  North  Dakota  Energy  &  Mineral  Research  Center 
noted  that  "domestic  production  currently  supplies  only  60%  of  our  petro¬ 
leum  requirements."  To  take  up  the  slack,  synthetic  fuels,  such  as  those 
derived  from  coal,  have  grown  increasingly  appealing  to  the  Department  of 
Defense  which,  according  to  Rindt,  "is  the  largest  single  consumer  of 
liquid  fuels  in  the  United  States."  The  by-products  of  coal  gasification 
are  a  potential  source  for  these  jet  fuels.  Therefore,  we  expect  that 
extensive  research  in  this  area  will  be  undertaken  to  define  and  improve 
the  thermal  stability  of  existing  jet  fuels  or  the  newly  coal-derived  jet 
fuels. 

When  jet  fuels  in  an  aircraft  system  are  exposed  to  elevated  tempera¬ 
tures,  it  causes  higher  thermal  stresses,  resulting  in  more  thermal 
degradation.  The  thermal  instability  manifests  itself  as  deposits  on  heat 
exchanger  surfaces,  fuel  control  valves  and  injector  feed  arms,  leading  to 
insufficient  heat  transfer  and/or  eventually  mechanical  problems.  Many 
experimental  studies^'^  have  been  performed  in  this  area.  No  single 
report  has  undertaken  a  detailed  coupled  thermal  fluids  modeling  effort, 
especially  a  3-D  analysis  of  flow  and  temperature  distributions,  for  the 
JFTOT^  experiments.  We  hope  this  report  will  provide  added  insight 
into  pertinent  parameters  that  need  to  be  treated  sel f-consi stently  to 
obtain  a  useful  computational  fluid  dynamics  and  chemistry  (CFDC)  model 
for  future  studies. 
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SECTION  III 


DEPOSITION  IN  LITERATURE 

This  section  contains  a  literature  survey*®  on  the  fouling  processes, 
chemical  kinetic  mechanisms  and  mathematical  models  of  deposition.  Fouling 
is  a  common  phenomena  occurring  in  many  industrial  processes:  chemical*7'2*, 
petroleum"4  aerospace43*,  nuclear43,  etc.  Fouling  that  consists  of  the 
deposit  of  corrosion  products  or  metal  ions  on  the  surface  is  classified  into 
the  corrosion  fouling*7,  particulate  fouling2®,  precipitation  fouling27, 
chemical  reaction  fouling  ,  etc.  Corrosion  fouling,  for  instance,  can 
manifest  itself  in  two  ways.  Corrosion  can  occur  at  the  heat-transfer  sur¬ 
face  and  the  resulting  corrosion  products  can  form  a  deposit  which  will 
inhibit  heat  transfer.  On  the  other  hand,  the  deposit  can  be  formed  by 
corrosion  products  produced  elsewhere  in  the  system,  which  are  then  trans¬ 
ported  to  the  heat-transfer  surface  and  deposited.  This  latter  form  of 
corrosion  fouling  is  an  occurrence  of  pariculate  fouling  or  of  precipita¬ 
tion  fouling.  The  processes  associated  with  given  categories  of  fouling 
appear  to  have  been  identified;  however,  in  general,  tneir  exact  character 
is  uncertain. 

In  nuclear  industries,  fouling  has  been  a  critical  issue  to  nuclear 
reactor  safety.2®  Deposit  of  particulate  metals  and  their  oxides  on  the 
cladding  material  of  the  nuclear  reactor  is  absolutely  detrimental  to 
reactor  safety  because  it  increases  core  temperatures  because  of  low 
thermal  conductivity  and  causes  significant  pressure  drops  in  a  narrow 
flow  channel  width,  leading  to  fuel-plate  instability.  In  order  to  over¬ 
come  this  problem  in  nuclear  reactors,  much  research  has  been  done  to 
define  the  deposit  mechanism  and  to  predict  the  deposit  thickness. 

Figure  1  shows  two  layers  of  oxide  deposited  on  the  aluminum  metal 
matrix  that  is  being  used  as  a  cladding  material  attached  to  uranium  oxide 
fuel.  The  deposition  is  basically  a  diffusion  process.  Depending  on  the 
interface  temperature,  aluminum  oxide  has  a  different  number  of  water 
molecules,  boehmite  ( A1 203  .  H20)  and  bayerite  (A1203.3H20) .  In  Figure  1, 
the  boehmite  layer,  the  inner  layer  on  the  surface,  grows  rapidly  to  the 
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with  Time 


Figure  1.  Cross  Section  of  Corroded  Aluminum. 
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maximum  thickness  while  the  outer  layer,  bayerite,  grows  in  thickness  with 
time,  and  as  the  flow  velocity  increases,  the  outer  layer  sloughs  off. 
Griess  et.al.2^  found  that  up  to  the  point  of  film  spallation  about  70% 
of  the  oxidized  aluminum  remained  on  the  surface  as  boehmite  and  the  rest 
was  lost  to  the  water. 

They  developed  an  empirical  deposit  correlation  considering  only  the  data 
acquired  with  the  water  at  a  coolant  pH  of  5.0  and  at  heat  fluxes  between 
3.15  and  6.3  x  106  W/m2: 

X  =  c  e0J78  exp  (-4600/T) 
where  X  =  oxide  thickness,  mils 

c  =  443  for  pH  5,  1200  for  pH  5.7  to  7,  constant 
6  =  exposure  time,  hours 

T  =  absolute  temperature,  oxide-water  interface,  degrees  Kelvin 


The  above  Arrhenius  type  equation  has  been  widely  used  to  predict  fouling 
deposits  on  the  heated  surface  of  nuclear  reactors  for  more  than  two 
decades. 

Regarding  the  chemical  kinetics  involved  in  deposition,  autoxidation 
kinetic  mechanisms  of  n-dodecane  are  well  defined  by  many  studies2^'2* 
to  explain  liquid  phase  oxidation  of  hydrocarbons  by  saturated  oxygen. 
Basically,  the  autoxidation  mechanism  is  as  follows: 

Initiation 


RH . Surface . >  R'  +  ’H 


(1) 
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Propagation 


R'  +O2 - >  RO2 ’ 


R02  •  +  RH - >  ROOH  +  R- 

Termination 


R02-  +  R02  *  - >  ROH  +  R1  -  CO  -  RJ  +  02 


R02-  +  R*  - >  ROOR 


R-  +  R-  - >  R  -  R 


(2) 

(3) 

(4) 

(5) 

(6) 


Initiation  occurs  on  the  hot  heated  surface  by  way  of  heterogeneous  cataly¬ 
sis  giving  alkyl  radicals  in  Reaction  1.  These  alkyl  radicals  react 
rapidly  with  molecular  oxygen  yielding  secondary  alkylperoxy  radicals, 

R02'  (Reaction  2).  Alkylperoxy  radicals,  on  the  other  hand,  react 
much  slowly  with  RH  to  form  al kyl hydroperoxides,  ROOH  shown  in  Reaction  3 
which  is  usually  the  rate  controlling  propagation  reaction.  Of  the  three 
termination  reactions,  Reaction  4  appears  to  be  the  dominant  process  since 
alcohols  and  ketones  are  detected  below  573  K.  Reaction  5  and  6  are  not 
important  because  of  low  concentration  of  R’.  The  autoxidation  reaction 
occurs  below  573  K.  Especially,  the  hydroperoxide  free  radical  in  step  3 
is  fairly  stable  under  560  K.  However,  most  experimental  data,  including 
Marteney  and  Spadaccini3^,  show  the  measurement  of  deposit  even  as  low 
as  500  K.  Also,  some  metal  ion  impurities  such  as  sulfur,  nitrogen,  etc. 
apparently  to  have  a  significant  effect  on  deposition.  Taylor33  carried 
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out  deposit  studies  using  six  different  jet  fuels  with  different  impuri¬ 
ties  and/or  with  deoxidation  treatment.  The  results  indicate  that  general 
trends  are  reduced  deposit  formations  with  deoxidized  fuels.  However,  a 
deoxidized  fuel  with  high  concentration  of  total  sulfur  gave  a  reversed 
result,  indicating  oxygen  is  not  the  only  element  that  causes  deposition. 
Also,  depending  on  the  chemical  structure  of  chemicals  associated  with 
sulfur  atoms,  sulfur  may  or  may  not  act  as  an  impurity.  For  instance, 
some  structures  of  sulfur  are  inert  even  in  the  deactivation  process  of 
chemical  catalysts.  Also,  the  analysis  of  deposits  indicates  that  the 
nature  of  deposits  formed  in  a  deoxygenated  fuel  is  different  from  those 
formed  in  an  air-saturated  fuel,  suggesting  that  such  deposits  are  formed 
as  a  result  of  different  classes  of  reactions.  Therefore,  detailed  reac¬ 
tions  leading  to  the  formation  of  actual  deposits  need  to  be  elucidated. 

Deposition  models  are  classified  into  reaction  controlled  deposi¬ 
tion34,  diffusion  controlled  deposition35,  and  a  combination  of  deposi¬ 
tion  and  removal36  by  shear  forces  or  turbulent  bursts. 

The  fouling  process  is  a  very  complicated  technical  area.  In  order  to 
model  the  fouling  process  completely,  mass  transfer,  fluid  dynamics,  and 
chemical  kinetics  in  models  are  required.  Studying  the  deposition  process 
for  oil -rigs  by  using  a  sand-water  slurry,  Watkinson  and  Epstein22  deve¬ 
loped  an  Arrhenius  type  correlation  to  determine  the  fouling  rate.  Using 
a  momentum  -  mass-transfer  analogy,  the  mass-transfer  coefficient  was 
related  to  fluid  velocity.  Also,  assuming  that  sticking  probability  of 
particles  on  the  surface  is  directly  proportional  to  adhesive  forces 
binding  a  particle  to  the  wall,  and  that  these  adhesive  forces  depend  on 
the  surface  temperature,  they  derived  the  following  Arrhenius-type 
equation. 


X 


A,(Cb  -  cw)  e'E/RT 
U«ff 


A2  f  U2  X 


(7) 
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where  Aj  and  A2  are  constants 

C  is  concentration  of  fouling  species 
U  is  velocity 

f  is  Fanning's  friction  factor. 


The  constants  Aj  and  A2  must  be  determined  for  a  specific  system 
and  operating  conditions.  Also,  one  thing  that  should  not  be  overlooked 
is  that  the  Schmidt  number  must  be  unity  in  order  to  make  a  momentum-mass 
transfer  analogy  valid.  Figure  2  contains  an  overall  view  of  fouling 
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process.  We  believe  that  chemical  reactions  are  taking  place  in  the 
bulk  stream  and/or  at  the  surface.  Products  generated  from  intermediate 
reactions  in  the  bulk  stream  move  toward  the  wall  by  diffusion. 

Somewhere  along  the  wall,  the  concentration  of  fouling  products  may 
become  greater  than  that  in  bulk  stream,  leading  to  a  back  diffusion. 
Also,  Figure  2  shows  the  deposit  removal  by  fluid  shear  or  turbulent 
bursts.  Generally,  diffusion  on  the  surface  is  defined  as  follows: 


Np  =  Kp  (Cpb  -  Cp.j ) 


or  assuming  the  first-order  chemical  reaction  on  the  surface 


Crittenden  and  Kol aczkowski^  have  extended  the  two-step  mass  transfer 
and  kinetics  model  to  include  back  diffusion.  Thus, 


dxf  =  .1 
dt  pf 


(Np  -  Nf) 


(8) 
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Figure  2.  Mass  Transfer,  Kinetics  and  Back  Convection  Model  Proposed  by 
Crittenden  and  Kolaczkowski 


where  Nf  =  Kf  (Cfi  -  Cfb) 


Assuming  »  Cfb  where  back  diffusion  occurs,  the  overall  fouling 
rate  excluding  removal  by  shear  is 


dRf  =  ! 

~5t  pf  kf 


(9) 


Replacing  mass-transfer  coefficients  and  reaction  velocity  constant  by 
fluid  velocities  and  the  Arrhenius  equation  respectively  gives  the  overall 
rate  as  follows: 


dRf  _  1  Cpb 

dt  Rfkf  p(d-2Xf)1,8  Sc0-67  i 

0.67XM0'2  G0-8  +  Aexp (-E/RT) 


0.607  \/i0-2G°'8Cf j 
P  (d~Xf )  1‘8SCf°-67 


M10) 


Equation  10  indicates  that  the  fouling  rate  depends  on  the  mass  flow  rate, 
G.  Also,  depending  on  the  system  and  operating  conditions,  it  can  be 
either  kinetics  or  diffusion  controlled,  leading  to  further  simplifica¬ 
tions  of  Equation  (10).  Including  deposit  removal  by  fluid  shear,  the 
overall  rate  can  be  written  as 


dRf 

~dtf  "  ”  ^2  Rf  t11) 


assuming  a  first  order  dependence  if  deposit  removal  rate  yields  an  asymp¬ 
totic  fouling  resistance-time  relationship  with  Rf  =  0,  when  t  =  0 
Rf  =  7r |  (1-  exp(-7T2  t)).  Kern  and  Seaton^®  proposed  that  deposit  removal  can 
be  assumed  to  be  proportional  to  wall  shear  stress  and  deposit  thickness: 
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(12) 


removal 


rate  = 


tfikf 


where  $  is  a  function  of  deposit  structure. 

As  described  above,  it  appears  that  many  deposition  models  are  avail¬ 
able  in  the  literature.  Parameters  in  those  models  are  system  dependent. 
No  information  appears  in  the  literature  about,  for  instance,  quantitative 
concentration  of  hydroperoxides  or  oxygen  in  any  jet  fuels  as  a  function 
of  temperature.  Therefore,  this  leads  us  to  use  a  single  correlation  on 
which  most  experiments  are  based. 
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SECTION  IV 


JFTOT  MODEL 

1.  NUMERICAL  SCHEME  AND  MODIFICATION 

Changes  to  the  KIVA  code  were  made  to  allow  for  the  representation  of 
fluids  other  than  an  ideal  gas.  These  changes  deal  primarily  with  the 
implicit  pressure  convergence  procedure  of  KIVA.^  The  partial  differen¬ 
tial  equations  solved  in  KIVA  are  the  Navier-Stokes,  conservation  of  mass, 
and  internal  energy  equations.  The  solution  procedure  for  these  equations 
in  KIVA  combines  the  Implicit  Continuous-fluid  Eulerian  ( ICE) ^  techni¬ 
que,  with  the  grid  rezoning  of  the  Arbitrary  Lagrangian-Eulerian  (ALE)^ 
method.  The  implicitness  of  this  method  is  the  iterative  updating  of 
portions  of  mass,  momentum,  and  internal  energy  conservation  equations, 
while  simultaneously  requiring  that  the  fluid  pressure  match  that  deter¬ 
mined  from  a  user  supplied  Equation  of  State  (EOS). 

KIVA  solves  the  fluid  equations  in  three  phases.  The  first  two, 
referred  to  as  Phases  A  and  B,  constitute  a  Lagrangian  calculation  in 
which  the  computational  grid  is  allowed  to  move  with  the  fluid.  During 
the  third  phase,  (Phase  C),  the  flow  field  variables  are  rezoned  onto  a 
new  grid.  Phase  A  is  an  explicit  updating  of  various  terms  of  the  fluid 
equations.  Phase  B  is  the  implicit  update  of  the  remaining  terms  of  these 
equations.  During  Phase  B,  the  following  finite-difference  approximations 
are  used: 

The  Pressure  Equation 


(13) 
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The  Volume  Equation 


vfjk  -  v?jk  +  At  I  (UA>; 


Equation  for  the  Cell -faced  Volume  Fluxes 


[«;  +  §;]  [(ua]= 


-  tuA]*  ]  =  -  At  I  [  [1  -  *,]p;  +  OqP5  +  I  P*X  ]  Aq .  Aa 


«[«' +  g«] »:  • 


Density  Equation 


v*.v 

A  '  Jk 


'  i  k  v  j  j  k 


Internal  Energy  Equation 


m?jk  [ink  -  i?ik]  =  - 


P?jk  +  P?Jk 


Kk  -  v*jk] 


The  superscript  of  these  equations  refers  to  the  value  of  that  particular 
variable  at  the  end  of  a  given  calculation  phase,  and  the  subscript  refers 
to  a  given  grid  location.  Some  of  the  finite-difference  approximations  of 
Reference  2  have  not  been  included  in  this  list,  since  they  do  not  enter 
into  the  present  discussion. 

In  order  to  converge  this  set  of  equations,  plus  those  not  listed 
above,  a  volume  residual  is  formed  from  Equation  (14)  as 


[  rv3jjk  =  vuk  -  [  v?,  +  At  z  (UA)  l  j 
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where  the  Phase  8  volume  is  derived  from  recasting  Equation  (13)  in  terms 
of  cell  pressure 


v 


B 

ijk 


(19) 


The  pressure  field  of  the  computational  grid  is  iteratively  updated  from 
Equations  15,  18,  and  19  by  the  Conjugate  Residual  algorithm  by  utilizing 
the  following  change  in  cell  volume  residual  with  respect  to  cell  pressure 


3rUk 


ap; 


ijk 


v  i  jk 

VA 

yi jk  i  Jk 


At2 

m  j  k 


(3*»  - 

4K  + 


(20) 


By  comparing  Equation  (20)  with  Equation  (18),  we  see  that  the  thermo¬ 
dynamic  volume  change,  because  of  pressure  for  this  iteration  scheme,  is 
as  follows 


av 

ap 


B 


ijk 

B 

ijk 


VA-. 

i  jk 


<yA  pA 
yi  jk*iJk 


(21) 


In  order  to  include  a  fluid  other  than  an  ideal  gas  in  this  conver¬ 
gence  scheme,  we  need  only  to  define  this  quantity  for  the  fluid  of 
interest.  For  the  present  application,  we  have  accomplished  this  by 
linearizing  the  EOS 


j  _  3P 

dP  "  a-p 


.dP  + 


dpi 

pT|p 


dl 


(22) 
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The  thermodynamic  volume  change  owing  to  pressure  can  be  defined  by 
differentiating  the  Density  Equation,  Equation  (16),  with  respect  to 
pressure. 


3V 

ap 


ills 

i  jk 


’tiXk  3pijk 

~7T~ 


i  jk 


V*-. 
v  i  j  k 


i  jk 


3Pijk 


i  jk 


(23) 


The  change  in  internal  energy  with  respect  to  pressure  during  this  phase 
of  the  calculation  can  be  obtained  by  differentiating  Equation  17  with 
respect  to  pressure,  with  the  added  approximation  that  the  average  pres¬ 
sure  of  the  P-V  work  term  equals  the  old  time  value.  The  final  result  is 
that  the  following  relationship  for  volume  changes  because  of  cell 
pressure 


3Vijk 


A  A  T  ?A 
M  j  k  V i i k  apl 


ap 


i  jk 


i  jk 


[ 


P; 


apli  /  1  ~  p2 


i  jk  3p 

31 


ijk 


] 


(24) 


The  only  remaining  task  is  to  define  the  EOS.  For  this  application 
the  following  Stiffened  EOS  was  adopted 


Pijk  as  (Pjjk  “  P0)  +  (Vijk  “  1  )  P  i  j  k 1 '  i  k 


(25) 


ap 

3P  I 


as  +  (>„k  "  ^.j. 


(26) 


ap 

ai 


Pi  jk  (yijk  ~  1  ) 

3s  +  (V„k  "  Ijl.J. 


(27) 


The  convergence  of  these  equations  was  tested  for  two  extreme  cases.  The 
first  was  that  of  an  ideal  gas.  The  results  were  essentially  identical  to 
those  of  the  original  KIVA  code.  The  second  was  in  the  limit  as  the  fluid 
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compressibility,  the  square  of  the  acoustic  velocity  of  the  stiffened  EOS, 
approaches  an  infinite  value  (that  is,  an  incompressible  fluid).  The 
results  of  this  test  indicate  that  as  the  velocity  was  increased,  the 
velocity  divergences  of  the  gird  cells  approached  zero,  as  would  be 
expected  for  an  incompressible  fluid. 

2.  GLOBAL  DEPOSIT  RATE 

The  composition  of  deposits  affords  clues  to  the  molecular  species 
involved  in  deposit  formation  and  the  mechanism  of  formation.  The 
chemical  reaction  mechanism  of  autoxidation  that  explains  deposition  can 
be  found  in  the  literature  as  mentioned  earlier.  However,  there  is  no 
single  datum  that  supports  the  quantitative  effect  of  those  chemical  free 
radicals  on  deposition  using  JFTOT. 

Since  the  majority  of  JFTOT  experiments  in  the  literature  were  not 
based  on  parameters  involved  in  chemical  reaction  or  mass  transfer,  a 
simple  global  type  deposit  rate  that  supports  most  JFTOT  experiments  was 
used  in  the  preliminary  calculation. 

The  experimental  data  set  used  in  obtaining  the  Arrhenius  factors  was 
from  Warner  and  Biddle.^  They  collected  real-time  deposition  rate  data 
using  the  Fiber  Optic  Modified  JFTOT  (FOM-JFTOT)  shown  in  Figure  3. 

Key  features  in  the  experiments  include  a  means  for  producing  quantita 
tive  data  with  the  JFTOT,  real-time  data  acquisition  and  display,  and  data 
at  up  to  eight  temperatures  obtainable  from  one  test.  Table  1  is  the  typi 
cal  data  set  of  deposits  that  was  built  up  to  0.14  microns.  Based  on 
Table  1,  the  Arrhenius  factors  were  obtained  using  a  linear  regression 
analysis. 


deposit  rate  =  155970  exp  ( - 1 9920/T ( K ) )  (28) 
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Figure  3.  Cross  Section  of  FOM-JFTOT  Heater  Tube  Holder. 


Table  1.  Real-Time  Deposition  Data  from  Warner  &  Biddle 


FF-3  JP-5  (470  B.P.) 

TEST  DATA  -  ALL  RUNS 

TIME  TO  ACHIEVE  1ST  MAX  (0.14u  DEPOSITS 


Test 

Tube 

Minutes 

Temp 

Probe 

Pos. 

Temp 

Run 

Run 

Run 

Run 

Run 

ID _ 

No. 

(mm) 

ILL 

No.  8 

No.  9 

No.  10 

No.  11 

No.  13 

550 

1 

28 

522 

107.5 

100.0 

92.5 

106.0 

2 

30 

531 

71.0 

63.0 

64.0 

74.5 

70.5 

3 

32 

538 

61.5 

52.0 

-- 

53.5 

58.5 

4 

34 

543 

44.0 

37.0 

46.5 

45.5 

47.5 

5 

36 

550 

38.0 

32.5 

37.5 

40.5 

38.5 

Run 

Run 

No.  15 

No.  16 

500 

1 

28 

475 

-  - 

2 

30 

482 

-- 

-- 

3 

32 

489 

245.0 

240.0 

4 

34 

495 

203.0 

-- 

5 

39 

500 

171.0 

156.0 

Run 

Run 

No.  27 

No.  28 

450 

1 

28 

426 

2 

30 

433 

-- 

-- 

3 

3? 

441 

-- 

4 

34 

444 

-- 

-- 

5 

39 

450 

763.3 

713.3 

18 


3.  JFTOT  RESULTS  ON  3-D  ANALYSIS 


Figure  4  shows  the  overall  grid  picture  for  the  JFTOT  model  of  8 
(radial)  x  24  (azimuthal)  x  24  (axial)  cells.  The  detailed  dimension  of 
Figure  4  is  in  Table  A-l.  All  the  figures  representing  the  JFTOT  model  in 
“us  report  are  not  to  exact  scale.  For  instance,  the  radial  and  azi¬ 
muthal  dimension  of  JFTOT  model  were  scaled  up  by  a  factor  of  5  to  provide 
a  clear  picture. 

Figure  5  is  the  top  view  of  8  (radial )x24(azimuthal )x24(axial )  cell 
model.  In  Figure  5,  the  azimuthal  nodes  of  the  inlet  and  outlet  tubes  to 
the  JFTOT  are  j=3  through  j=7  for  inlet,  and  j=9  through  j  =  1 3  for  outlet. 
The  90°  rotation  and  the  cross-sectional  area  of  the  inlet  and  outlet 
tube  were  simulated  to  be  identical  to  the  JFTOT  dimension  (4.901e-6  m2, 

5  azimuthal  nodes,  and  2  axial  nodes  in  model). 

The  inner  boundary  is  the  heater  surface  of  the  JFTOT.  The  surface 
temperature  distributions  on  the  heater  were  obtained  from  Warner  & 
Biddle's  Table  (550  F  case).  Since  Warner  &  Biddle's  experiments  included 
only  part  of  the  test  section  (2.8  cm  to  3.9  cm),  the  predicted  tempera¬ 
ture  profile  was  correlated  using  a  regression  analysis.  The  following 
correlation  was  incorporated  into  KIVA. 

T(K)  =  519.6  -  3.676X  +  0.2718X2  -  0.003857X3  (29) 

where  X  is  axial  length  in  mm. 

The  outer  boundary  temperature  was  taken  as  300  K  for  the  entire  length  of 
JFTOT.  The  inlet  velocity  to  the  JFTOT  was  set  at  0.01  m/s. 

Figure  6  (2600  cycles,  time  =  4.1  sec)  shows  the  temperature  contours 
and  velocity  vectors  for  the  three  azimuthal  slices  corresponding  to  the 
outlet  plane,  mid-height,  and  exit' plane.  The  temperature  contours  show 
the  penetration  of  thermal  energy  into  the  flow.  At  the  inlet  location, 
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Figure  4 


Overall  Grid  View  of  JFTOT  Model  of  8x24x24  Cells 
(Base  Case) . 
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Figure  5.  Top  View  of  JFTOT  Model  of  8x24x24  Cells. 
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Figure  6.  Calculated  Temperature  Contours  and  Velocity  Vectors 
at  the  Bottom,  Mid-Height,  and  Top  of  the  FOM-JFTOT 
(2600  Cycles,  Time  =  4.1  Sec). 


the  temperature  contours  are  closely  spaced  in  the  vicinity  of  the  heater 
tube  surface  on  the  side  of  the  tube  closest  to  the  inlet.  On  the  oppo¬ 
site  side  of  the  tube  the  contours  are  more  widely  spread. 

This  indicates  that  thermal  energy  is  penetrating  into  the  fuel  much 
further  on  the  back  side  of  the  tube,  where  radial  and  azimuthal  veloci¬ 
ties  are  nearly  stagnant.  At  the  mid-height  location,  contours  are  almost 
evenly  spaced,  while  the  velocity  vectors  appear  to  be  outward.  However, 
since  vectors  are  originated  inside  the  computing  grid  domains,  this  indi¬ 
cates  the  instantaneous  direction  of  the  local  flow.  A  better  picture  can 
be  seen  in  the  Appendix.  The  maximum  axial  and  azimuthal  velocity  indexed 
as  urn  and  vm,  respectively,  are  much  smaller  than  that  of  axial  velocity 
(wm)  by  more  than  a  factor  of  200. 

As  we  move  to  the  azimuthal  plane  at  the  exit,  smaller  temperature 
gradients  are  apparent  on  the  side  corresponding  to  the  flow  exit,  while 
in  the  other  region,  velocity  contours  are  nearly  symmetric.  The  magni¬ 
tude  of  the  velocity  vectors  are  not  relatively  comparable  from  one  plane 
to  the  others. 

KIVA  calculates  fluid  temperatures  corresponding  to  each  computational 
cells  from  minimum  to  maximum  indexed  as  min  and  max,  respectively  in 
Figure  6.  However,  temperature  contours  only  between  1  (low)  and  h  (high) 
are  plotted  in  this  figure. 

Figure  7  shows  the  radial  profiles  of  temperature  and  velocity  at  the 
same  axial  plane  as  Figure  6  for  the  azimuthal  location  corresponding  to 
the  exit.  The  velocity  profiles  at  the  entrance  and  mid-height  are  basi¬ 
cally  parabolic,  indicating  a  simple  pipe  flow,  whereas  at  the  exit  the 
axial  velocity  at  the  outer  boundary  indicates  the  presence  of  the  outlet 
pipe.  The  axial,  radial,  and  azimuthal  velocities  at  the  exit  are  nearly 
in  the  same  order  of  magnitude  as  shown  in  Figure  6.  At  the  inlet  plane, 
the  temperature  gradients  appear  steeper  than  those  of  mid-height,  and 
exit.  A  steep  temperature  gradient  near  the  tube  surface  is  equivalent  to 
the  dense  contours  shown  in  Figure  6.  At  the  mid-height  the  temperature 
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Radial  Position  (Nodes)  Radial  Position  (Nodes) 


Figure  7.  Calculated  Radial  Temperature  and  Velocity 
Distributions  at  the  Bottom,  Mid-Height, 
and  Top  of  the  FOM-JFTOT. 


Axial  Velocity,  m/s 


gradient  is  still  steep,  but  not  as  steep  as  it  was  at  the  inlet  plane. 
This  again  corresponds  to  the  spreading  of  the  temperature  contour  shown 
in  Figure  6.  At  the  exit  plane,  the  temperature  gradient  is  not  signifi¬ 
cantly  large  across  the  tube  radius,  indicating  that  the  thermal  energy 
has  penetrated  the  flow  channel  gap. 

An  added  view  of  the  flow  is  presented  in  Figure  8,  which  shows  the 
entire  axial  velocity  vectors  corresponding  to  the  azimuthal  slice  of 
inlet  and  exit  planes,  respectively. 

At  the  inlet  plane,  the  velocity  vectors  on  the  opposite  side  of  the 
inlet  stagnate.  At  the  mid-height,  the  velocity  vectors  appear  symme¬ 
tric.  The  complex  radial  and  azimuthal  mid-plane  flow  pattern  shown  in 
Figure  6  was  not  significant  to  these  calculations  because  their  magni¬ 
tudes  are  much  smaller  than  those  of  Figure  8.  At  the  exit,  vectors  on 
the  opposite  side  of  the  exit  appear  small  as  would  be  expected. 

3.1  Sensitivity  Study  of  Computational  Grids 

The  following  two  cases  were  examined  to  find  a  sensitivity  of  added 
or  reduced  computational  nodes  on  the  predicted  temperature  and  velocity 
distributions  described  in  the  preceding  section,  hereafter  called  the 
Base  Case. 

The  computational  cells  or  elements  of  Case  1  consist  of  8(radial)x 
24(azimuthal ) x 1 4 ( ax i al ) ,  which  is  a  reduced  number  of  cells  in  the  axial 
coordinate  from  that  of  the  Base  Case.  All  cases  for  the  sensitivity 
studies  were  run  up  to  2600  cycles,  time  =4.1  seconds. 

Figure  9  shows  the  overall  grid  for  Case  1.  Again,  the  dimension  of 
the  radial  and  azimuthal  coordinates  was  scaled  by  a  factor  of  5.  In 
orde"  to  have  the  same  cross-sectional  area  of  the  inlet  and  the  outlet 
tube  as  that  of  the  base  case  (4.91  e-6  m2),  the  height  of  the  inlet  and 
outlet  tube  remained  the  same  as  that  of  the  base  case.  The  detailed 
dimension  «f  the  computational  grids  is  given  in  the  Appendix. 
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Figure  10  is  equivalent  to  Figure  6.  The  temperature  contour  pattern 
appears  nearly  the  same  as  those  in  Figure  6.  Since  the  surface  tempera¬ 
ture  correlation  (eq.  29}  is  a  parabolic  profile,  the  less  accurate  axial 
description  resulted  in  a  lower  fluid  exit  temperature.  The  velocity  vec¬ 
tors  at  the  mid-height  appear  inward  especially  on  the  side  corresponding 
to  the  inlet  tube,  while  the  velocity  vectors  in  Figure  6  are  more  dis¬ 
persed  at  the  mid-point.  This  indicates  that  with  more  computational 
cells  in  the  axial  coordinates,  radial  and  azimuthal  flow  patterns  appear 
more  dispersed.  However,  the  fairly  insignificant  magnitude  of  the  radial 
and  azimuthal  velocity  compared  to  that  of  the  axial  velocity  does  not 
change  the  thermal  equilibrium. 

Figure  11  contains  the  radial  temperature  and  axial  velocity  profiles 
at  three  axial  locations,  equivalent  to  those  of  Figure  7.  No  differences 
are  noted  between  those  two  figures. 

Figure  12  shows  the  entire  axial  velocity  vector  profile  azimuthally 
sliced  at  the  location  corresponding  to  the  inlet  and  outlet  tube,  which 
can  be  compared  to  Figure  8.  Some  vectors  at  the  mid-height  shown  in 
Figure  12  show  inward  flow  pattern  as  also  shown  in  Figure  10.  However, 
as  mentioned  earlier,  the  magnitude  of  the  radial  and  azimuthal  velocity 
vectors  is  much  smaller  compared  to  those  of  axial  velocity  to  affect  the 
thermal  energy  distribution  in  JFTOT. 

The  computational  cells  of  Case  2  shown  in  Figures  13  and  14  consist 
of  8(radial )xl2(azimuthal )xl4(axial )  cells,  which  is  a  reduction  in  the 
number  of  cells  in  the  azimuthal  direction  compared  to  Case  1.  Again,  the 
cross-  sectional  area  of  the  computational  cells  was  identical  to  the  size 
of  inlet  and  outlet  tubes  of  F0M-JFT0T.  Figure  15  shows  the  obtained 
temperature  contours  and  velocity  vectors  azimuthally  sliced  at  the  inlet, 
mid-height,  and  exit.  These  slices  compare  with  Figures  6  and  10  for  the 
Base  Case  and  Case  1,  respectively. 
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Radial  Position  (Nodes) 


Figure  11.  Calculated  Radial  Temperature  and  Velocity 
Distributions  at  the  Bottom,  Mid-Height,  and 
Top  of  the  FOM-JFTOT  (Case  1). 
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Figure  15.  Calculated  Temperature  Contours  and  Velocity 

Vectors  at  the  Bottom,  Mid-Height,  and  Top  of  the 
FOM-JFTOT  (Case  2,  2600  Cycles,  Time  =  4.1  Sec). 


We  see  no  significant  differences  among  those  three  cases  except  the 
velocity  vector  profile  at  the  mid-height.  With  more  computational  cells 
either  in  the  azimuthal  or  in  the  axial  coordinate,  the  radial  and 
azimuthal  velocities  appear  more  dispersed  or  damped.  For  instance,  the 
ratios  of  the  axial  velocity  to  the  radial  velocity  at  the  mid-height  are 
241,  33,  and  25  for  the  Base  Case,  Case  1,  and  Case  2,  respectively. 
However,  the  thermal  energy  distributions  are  dominated  by  axial  flow 
motions. 

Figure  16  is  comparable  to  Figures  7  and  11.  Radial  temperature  and 
axial  velocity  profiles  corresponding  to  the  exit  are  nearly  identical  to 
those  of  Figures  7  and  11. 

Figure  17  is  the  axial  velocity  profile  azimuthally  sliced  at  the 
inlet  and  the  opposite  side  of  the  inlet,  and  the  exit  and  the  opposite 
side  of  the  exit  that  can  be  compared  to  Figures  8  and  12.  Figure  17 
shows  that  the  axial  velocity  at  the  inlet  and  exit  for  the  reduced 
computational  grid  is  nearly  the  same  as  those  for  a  more  refined  grid, 
indicating  that  the  computational  grid  of  8(radial )xl2(azimuthal )x 
1 2 ( axi al )  is  fairly  good  at  simulating  the  JFTOT,  thereby  avoiding  a  lot 
of  computer  costs. 

3.2  Sensitivity  Study  of  Computing  Time 

The  base  case  of  8  x  24  x  24  cells  was  run  up  to  3950  cycles  to  see  if 
the  temperature  and  velocity  had  reached  a  steady  state.  Figures  15  and 
16  were  velocity  profiles  of  3650  cycles,  time  =  5.8  seconds  and  3950 
cycles,  time  =  6.42  seconds  at  the  mid-height.  These  are  compared  with 
the  velocity  vector  profile  at  the  mid-height  shown  in  Figure  6. 

As  noted  in  Figures  6,  18,  and  19,  velocity  vectors  are  more  azimu¬ 
thally  oriented  as  time  goes  on.  Also  one  can  note  an  azimuthally  oscil¬ 
lated  flow  pattern  at  the  different  computational  time  frame,  indicating 
that  azimuthal  and  radial  velocities  are  still  transient.  However,  the 
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Figure  16.  Calculated  Radial  Temperature  and  Velocity  Distributions 
at  the  Bottom,  Mid-Height,  and  Top  of  the  FOM-JFTOT 
(Case  2). 


Velocity,  m/s 


The  Axial  Velocity  Vectors  Corresponding  to  the 
Azimuthal  Slices  of  the  Inlet  and  Exit  Plane  (Case  2) 


um=  9.97401  e-03  vm=  9.1 8760e-03  wm=  4.60424e-01 


18.  Calculated  Velocity  Vectors  at  the  Mid-Height 
for  Sensitivity  Study  (Base  Case,  3650  Cycles, 
Time  =5.8  Sec) . 


um=  7.58489e-03  vm=  5.43037e-03  wm=  3.9931 8o-03 


Figure  19.  Calculated  Velocity  Vectors  at  the  Mid-Height  for 
Sensitivity  Study  (Base  Case,  3950  Cycles, 

Time  =  6.42  Sec) . 


temperature  contours  are  dominated  by  the  axial  velocity  and  Figures  20 
and  21  indicate  that  temperatures  appear  to  be  at  steady  state. 

3.3  Deposit  Calculation 

The  inner  aluminum  tube  surface  temperature  (eq.  29)  was  incorporated 
into  KIVA,  along  with  the  physical  property  correlation  of  JP-5.  The 
azimuthally  uniform  temperature  and  the  interface  temperature  300  K 
between  the  fluid  and  the  inner  housing  were  assumed.  The  system  pressure 
and  volumetric  flow  rate  used  in  this  computation  were  3.45  MPa  and 
0.003  liter/min,  respectively. 

The  rate  of  formation  of  deposit,  dx/dt,  is  defined  by  the  following 
equation  as  mentioned  earlier  in  the  Section  IV  (subsection  2). 


dx/dt  =  155970  e  (-19920/T),  m/sec 


where  T  is  the  temperature  of  the  tube  wall  expressed  in  degrees  Kelvin. 

Boundary  heat  flux  and  fluid  temperatures  obtained  from  KIVA  were  used 
as  boundary  conditions  to  calculate  the  azimuthal  wall  temperature  and 
carbon  deposit  thickness  by  solving  the  following  energy  equation  sepa¬ 
rately  from  KIVA  in  conjunction  with  the  carbon  deposition  correlation 
above. 


KV2T 


0 


(30) 


Implicit  in  a  deposition  model  such  as  this  are  a  number  of  assump¬ 
tions  about  the  deposition  process.  First,  it  is  assumed  that  the  deposi¬ 
tion  depends  on  the  temperature  at  the  heater  wall.  It  is  also  assumed 
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nLn=  3. 1 1957E+02  mx=  5.06272^02  1=  3.31380E+O2  h=  4.86840E*D2 


Figure  20.  Calculated  Temperature  Contours  at  the 
Mid-Height  for  Sensitivity  Study 
(Base  Case,  3650  Cycles,  Time  =  5.8  Sec). 


3. 12271E+02  m a 


1=  3.31901E+02  h=  4.889460  02 


'e  21.  Calculated  Temperature  Contours  at  the  Mid-Height 
for  Sensitivity  Study  (Base  Case,  3950  Cycles, 

Time  =  6.42  Sec). 
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that  the  deposits  accumulated  at  any  position  along  a  wall  are  due  to  the 
local  wall  temperature  and  are  not  affected  by  the  transport  mechanism  in 
the  vicinity  of  that  particular  location.  Data  to  model  the  formation  and 
transport  of  particles  in  the  bulk  do  not  exist,  so  global  models  have  to 
suffice  for  initial  studies. 

Figure  22  shows  axial  profiles  of  the  wall  temperature  and  deposit 
thickness  for  the  azimuthal  location  corresponding  to  the  exit.  The  model 
being  employed  indicates  that  the  growth  of  the  deposit  layer  clearly 
follows  the  growth  of  the  temperature  along  the  tube.  The  maximum  deposit 
thickness  of  0.13  microns  was  calculated  at  38  min  at  the  hot  spot  where 
the  wall  temperature  is  560  K,  which  compares  favorably  with  the  maximum 
measured  deposit  of  0.14  microns.  Furthermore,  the  model  is  consistent 
with  the  experiment  in  predicting  the  location  of  the  maximum  deposit  at 
the  location  of  maximum  tube  temperature.  However,  the  validity  of  this 
model  is  based  on  the  very  limited  existing  experimental  data  and  assump¬ 
tions  made  for  the  JFTOT  boundary  condition.  This  is  primarily  due  to  a 
lack  of  data  of  this  nature.  However,  it  is  a  start. 
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SECTION  V 


CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  STUDIES 


Temperature  and  velocity  distributions  in  the  FOM-JFTOT  calculated 
from  this  study  offer  an  aid  in  understanding  fluid  behavior  that  would 
eventually  need  to  be  coupled  with  the  chemistry  and  heat  transfer  later 
in  a  future  study.  The  CFDC  model  used  in  this  study  is  based  on  the 
global  Arrhenius  type  surface  temperature-dependent  correlation  because  of 
the  lack  of  the  existing  experimental  data.  The  initial  attempt  to 
simulate  flow  in  the  JFTOT  is  very  encouraging.  "Clean"  experiments  need 
to  be  done  to  define  all  the  parameters  pertinent  to  the  deposition.  The 
following  are  recommended  for  future  studies. 

1.  Quantitative  measurement  of  free  oxygen  and  hydroperoxide  radicals  in 
the  bulk  stream  as  a  function  of  temperature. 

2.  Accurate  flow  distribution  measurement. 

3.  Measurement  of  JFTOT  housing  temperature  (interface  temperature 
between  fluids  and  housing). 

4.  Qualitative  data  acquisition  on  particle  removal  from  the  deposit 
layer  as  a  function  of  temperature  and  flow  rate. 

5.  Qualitative  (for  instance,  what  chemical  structure  of  sulfur)  and 
quantitative  impurity  effect  on  deposition. 
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NOMENCLATURE 


speed  of  sound  in  fluid 
constant  defined  in  eqn.  7 
constant  defined  in  eqn.  7 
outward  projected  area 
concentration  of  fouling  species 
pH  constant 

activation  energy  defined  in  equation  7 

friction  factor 

mass  flux 

internal  energy 

reaction  rate  constant 

mass  transfer  coefficient 

thermal  conductivity  defined  in  eqn.  30 

cell  face  mass 

mass  flux 

pressure 

volume  residual 

fouling  resistance 

face-centered  coupling  coefficient 

Schmidt  number 

time 

temperature 

velocity 

velocity 

volume 

deposit  thickness 

axial  distance  defined  in  equation  29 
deposit  rate  defined  in  eqn.  7 


face  number 

exposure  time 

compressibility 

defined  from  f  =  ARe“0-2 

viscosity 

defined  in  eqn.  10 
defined  in  eqn.  11 
density 

wall  shear  stress 

variable  implicitness  parameter 

function  of  deposit  structure 
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SUPERSCRIPT 

A  phase 

B  phase 

SUBSCRIPT 

b  bulk 

B  bottom  faces  of  cell  ijk 

D  derriere  faces  of  cell  ijk 

f  foul  ant 

f^  foul  ant  at  bulk 

fi  foulant  at  interface 

F  front  faces  of  cell  ijk 

i  fluid-deposit  interface 

L  left  faces  of  cell  ijk 

m  mixture 

o  referenced  state 

p  precursor 

Pb  precursor  at  bulk 

pi  precursor  at  interface 

R  right  faces  of  cell  ijk 

T  top  faces  of  cell  ijk 

w  wall 
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APPENDIX  A 


DIMENSION  TABLES  FOR  JFTOT  MODEL 


Table  A-l.  Axial  Dimension  of  8x24x24  Cells 


Axial  Distance  (ml 

K  = 

1 

0.0000e-3 

K  = 

2 

1 . 0138e-3 

K  = 

3 

2 . 0276e-3 

K  = 

4 

4.8249e-3 

K  = 

5 

7.6221e-3 

K  = 

6 

1.04194e-2 

K  = 

7 

1.32166e-2 

K  = 

8 

1 . 601 39e - 2 

K  = 

9 

1 . 881 1 le-2 

K  = 

10 

2 . 16084e-2 

K  = 

11 

2.44056e-2 

K  = 

12 

2 . 72028e - 2 

K  = 

13 

3 .00000e-2 

K  = 

14 

3.27973e-2 

K  = 

15 

3 . 55945e - 2 

K  = 

16 

3.83918e-2 

K  = 

17 

4.11 890e - 2 

K  = 

18 

4. 39863e-2 

K  = 

19 

4.67835e-2 

K  = 

20 

4 . 95808e-2 

K  = 

21 

5.23780e-2 

K  = 

22 

5. 51753e-2 

K  = 

23 

5. 79725e-2 

K  = 

24 

5.89862e-2 

K  = 

25 

6 . 00000e-2 

A- 1 


Table  A-2.  Axial  Dimension  of  8x24x14  Cells 


Axial  Distance  (ml 

K  = 

1 

0.00Q0e-3 

K  = 

2 

1.0138e-3 

K  = 

3 

2.0276e-3 

K  = 

4 

7.6221e-3 

K  = 

5 

1.32166e-2 

K  = 

6 

1 . 881 1 le~2 

K  = 

7 

2.44056e-2 

K  = 

8 

2.00000e-2 

K  = 

9 

3 . 55945e-2 

K  = 

10 

4. 1 1890e - 2 

K  = 

11 

4 . 6783  5e - 2 

K  = 

12 

5.23780e-2 

K  = 

13 

4.79725e-2 

K  = 

14 

5.89862e-2 

K  = 

15 

6.00000e-2 

K  =  1,  2,  and  3  open  for  inlet 
K  =  13,  14,  and  15  open  for  exit. 


Table  A-3.  Radial  Dimension  of  JFTOT  Model 


i  =  1 
i  =  2 
i  =  3 
i  =  4 
i  =  5 
i  =  6 
i  =  7 
i  =  8 
i  =  9 


Radial  Distance  (m) 

1.4000e-3 
1 . 5625e-3 
1 . 6696e - 3 
1 . 7768e-3 
1.8839e-3 
1 . 99 1 1 e - 3 
2.0982e-3 
2.2054e-3 
2.3125e-3 


i  =  1  is  the  centerline  of  the  heater  rod  thickness  for  the  inner 
boundary. 

Computational  grid  for  fluid  is  from  i  =  2  to  i  =  9. 
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Table  A-4.  Azimuthal  Dimension  of  JFTOT  Model 


For  12  cells,  0.5236  radians 
j  =  2,  3,  and  4  for  inlet 
j  =  5,  6,  and  7  for  exit. 

For  24  cells,  0.2618  radians 
j  =  3,  4,  5,  6,  and  7  for  inlet 
j  =  9,  10,  11,  12,  and  13  for  exit. 
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program  cfilm 

This  is  carbon  thickness  calculation  for  the  final  report. 


c 
c 
c 

c  the  unit  is  in  S.I. 

c  heat  flux  =  w/sq.m,  thermal  conductivity  =  w/m-K  . 
c  temperature  =  K,  heat  transfer  coefficient  =  w/sq.m-K  .  c 

tkiva  =  fluid  bulk  temperature  next  to  the  wall,  K. 
c  hkiva  =  convective  heat  transfer  coefficient,  w/m*m-K.  c 

condcl  =  thermal  conductivity  of  carbon  deposit,  w/m-K.  c 
condhr  =  thermal  conductivity  of  the  wall,  w/m-K. 
c  qkiva  =  average  heat  flux  calculated  from  KIVA,  w/m*m.  c 

aeon  =  pre-exponential  constant, 
c  bcon  =  exponential  constant, 

implicit  real*8 (a-h, o-z) 

dimension  t (5, 24) , condi (5,24) , condj (5,24) , aw (5, 24) ,ae(5,24) 

. ,  an (5 , 24) , as ( 5 , 24 ) , ap ( 5 , 24) , hkiva (24) , qbc (24 ) 

.,  rad (5) , tkiva (24) , r (5, 24) , fthk(24) ,qrad(5, 24)  data 

condcl  /  0.121/ 

data  condhr  /  18.83/ 

data  omega  /  1.5/ 
c 

c  k  =  15 

c 

data  tkiva  /  506.,  506.,  506.,  506.,  506.,  506., 

+  506.,  506.,  506.,  506.,  506.,  506., 

+  507.,  507.,  507.,  507.,  507.,  507., 

+  507.,  507.,  506.,  506.,  506.,  506.  / 

data  hkiva  /  12852.,  12852.,  12852.,  12852.,  12852.,  12852., 

+  12852.,  12852.,  12852.,  12852.,  12852.,  12852., 

+  13094.,  13094.,  13094.,  33094.,  13094.,  13094., 

+  13094.,  13094.,  12852.,  12852.,  12852.,  12852. 

/  data  qkiva  /  6.94e5/ 

data  aeon  /  1.5597e5/ 
data  bcon  /  19920./ 
data  timend  /  2250./ 
data  dt  /  60.0  / 

data  rad  /  0.001150,  0.00135625,  0 . 0015625 , 0 ., 0 . / 
c 

open (6 , f ile= ' cfilm. out ' , status=  *  unknown ' ) 
dthko  =  l.e-6 
rad (4)  =  rad (3)  +  dthko 
rad (5)  =  rad (4)  +  dthko 
delth  =  2. *3. 14159/24 
qkiva  =  qkiva*rad(5)/rad(l) 
c 

c  *  *  *  form  conductance  matrix 
c 

do  10  i  =  1,5 

do  10  j  =  1,24 

r(i, j )  =  rad (i) 

condi (i,j)  =  condhr 

if(i.eq.4)  condi(i,j)  =  condcl 

if (i.eq.5)  condi(i,j)  =  hkiva (j )* (rad (i)  -  rad(i-l)) 
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10  continue 
c 

do  20  i  =  1,5 
do  20  j  =  1,24 
cond j ( i , j )  =  condhr 
if(i  .eq.  4)  condj(i,j)  =  condcl 
if(i.eq.5)  conej(i,j)  =  condcl 
20  continue 
c 

do  30  i  =  1,5 
do  30  j  =  1,24 
t(i,j)  =  tkiva(j) 

30  continue 
c 

100  time  =  time  +  dt 
c 

do  110  j  =  1,24 
c 

c  temp  =  561  K  case 
c 

r(4,j)  =  r (4 , j )  +  dt*acon*exp(-bcon/t (4, j) ) 
r(5, j)  =  r (4 , j )  +  dthko 
110  continue 
c 

do  120  i  =  1,4 
do  120  j  =  1,24 
im  =  i-1 

if (im.eq. 0)  im  =  1 
ip  =  i  +  1 

rp  =  0.5*(r(ip,j)+r(i,j)) 

rm  =  0.5*(r(im, j)+r(i, j) ) 

drp  =  amaxl ( (r ( ip , j )  -  r(i, j) ), dthko) 

drm  =  amaxl ( (r(i, j)  -  r(im, j) ) , dthko) 

dr  =  rp  -  rm 

an(i,j)  =  condj (i, j)/(r(i, j) *delth) **2 

as(i,j)  =  an(i, j ) 

ae(i,j)  =  condi ( ip , j ) *rp/ ( r ( i , j ) *drp*dr ) 
if(i.eq.l)  then 
aw(i,j)  =  o.O 

qbc(j)  =  rm*qkiva/(dr*r(i, j) ) 
else 

aw(i, j)  =  condi ( i , j ) *rm/ ( r ( i , j ) *drm*dr ) 
endif 

ap(i,j)  =  aw(i, j)+ae(i,j)+as(i, j)+an(i, j) 
120  continue 
c 

nloop  =  o 

1000  nloop  =  nloop  +  l 
deltmx  =  0.0 
c 

do  130  i  =  1,4 
do  130  j  =  1,24 
c 

im  =  i-i 
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if(i.eq.l)  im  =  i 
ip  =  i+1 
jP  =  j+1 

if (j  *eq.24)  jp  =  1 
jm  =  j-1 

if ( j . eq. 1  )  jm  =  24 
fig  =  0.0 

if(i.eq.l)  fig  =  l.o 
tsave  =  t(i,j) 

t(i,j)  =  (aw(i, j) *t(im, j)+ae(i, j)*t(ip, j)+an(i, j)*t(i,  jp) 
+as(i,] ) *t (i, jm) +f lg*qbc ( j ) )/ap (i ,  j ) 
t(i, j)  =  tsave  +  omega* (t(i,j)  -  tsave) 
delt  =  abs(t(i,j)  -  tsave) 
if (t(i, j) .ne.0.0)  delt  =  de!lt/t(i,j) 
deltmx  =  amaxl (deltmx, delt) 

130  continue 

if (deltmx. gt. 5. e-6  .and.  nloop. It . 20000)  go  to  1000 

if (amod(time, 60. *10. )  .ne.  0.  .and. 

•  time  .It.  timend  )  go  to  160 

do  135  i  =  1,4 
do  135  j  =  1,24 
ip  =  i+1 

qrad ( i , j )  =  -condi (ip, j) * (t (ip, j) -t(i, j ) )/ (r (ip, j ) -r (i, j ) ) 

135  continue 


timem  =  time/ (60.) 

write (6,*)  'time (min) , timem 

do  140  j  =  1,24 

write  (6,  *)  'j=\j 

write (6, 146)  (t(i, j) ,  i=l,5) 

write (6, 147)  (qrad (i, j ) ,  i=l,4) 

140  continue 

145  format(5(lx,lpell.4) ) 

146  formate  T  =' ,5(lx,  lpell.4) ) 

147  formate  qr= »  ,  5  (lx,  lpell .  4 ) ) 
do  150  j  =  1,24 

fthk(j)  =  r(4, j)  -  r(3, j)  -  dthko 
150  continue 

write (6,*)  '  carbon  film  thickness' 
write (6, 145)  (fthk(j),  j=l,24) 

160  if (time. It. timend)  go  to  100 

stop 

end 
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APPENDIX  C 
ENLARGED  FIGURES 
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